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T-matrix  Approach  to  Array  Modeling 


by 

C.  L.  Scandrett  (Department  of  Mathematics) 
S.  R.  Baker  (Department  of  Physics) 
Naval  Postgraduate  School 
Monterey,  CA 

Abstract 


The  so-called  T-matrix  formulation  for  transducer  array  calculations  is  described  from 
a  theoretical  standpoint.  The  methodology  is  then  applied  to  an  array  problem  whose 
solution  can  be  found  analytically,  thereby  producing  a  reference  solution  for  comparing 
array  calculations  based  upon  FEM/BEM  determined  T-matrices.  Three  numerical  codes 
have  been  tested  for  their  accuracy  in  reproducing  the  exact  T-matrix  for  a  spherical  shell. 
These  in  turn,  are  used  to  calculate  the  near  and  far  field  pressure  for  the  reference  array. 
It  is  demonstrated  in  the  report  that  the  methodology  appears  to  be  robust.  Provided  an 
“accurate”  finite/boundary  element  model  of  a  transducer  is  used  to  produce  the  transducer’s 
T-matrix,  “accurate”  field  pressure  results  for  both  the  far  and  near  field  of  the  array  can 
be  found. 
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1  Theoretical  Background 

The  time  harmonic  (e!wi)  pressure  field  radiated  or  scattered  from  a  fluid  loaded  finite  struc¬ 
ture  can  be  represented  in  terms  of  outgoing  spherical  Hankel  functions  applied  to  the  set  of 
spherical  harmonics.  For  a  single  radiator/scatterer  with  a  local  coordinate  system  written 
with  the  index  “j”,  the  functional  form  of  the  pressure  field  (pj)  can  be  written 

N  n 

Pj(rj,0j,<l>j)  =  Y  Aimnhn{krj)Sl™(eh<t>j) 

n= 0  m~—n 

where 

=  P?(cosO)eim* 

are  spherical  harmonic  functions. 

A  distinction  is  made  between  radiated  pressures  of  a  structure  Bjmn  and  scattered 
pressures  Ajmn  throughout  the  analysis  used.  It  is  assumed  that  the  coefficients  for  the 
radiated  amplitudes  are  known  or  have  been  calculated  using  the  finite  element  method, 
while  the  scattered  pressures  are  to  be  determined  when  a  series  of  transducers  are  placed 
in  an  array,  and  are  insonified  by  either  an  exterior  source  or  from  mutual  array  element 
interactions  when  the  array  is  in  an  active  mode  of  operation. 

What  the  spherical  addition  theorem  allows  one  to  do,  is  represent  a  pressure  field  relative 
to  one  coordinate  system  in  terms  of  coordinates  of  a  second  system.  The  representation  of 
a  single  outgoing  spherical  wave  is  given  by  the  formula  [3] 

oo  v  n+is  f 

K(kr ,)«?(«.,*)  =  EE  £ 

„=0  ti--v  p=|n_t,| 

p>\m-n\ 

a(v,p,n,p,,  m)jl/(kr2<)hp(kr2>  )Q£+TO  (02l ,  <f>21  (02,  <f>2 ) 

(the  prime  on  the  summation  over  p  indicates  jumps  of  2  in  the  sum) 

(ri,6i,(f>i)  —  spherical  coordinates  relative  to  system  1 
(r2,62,<f>2)  =  spherical  coordinates  relative  to  system  2 
(r2i,  @2\ 1 4>2\)  —  origin  of  system  2  relative  to  1 
r2>  =  max{r2,r2i}  r2<  =  min{r2,r2i} 

and  where  the  coefficients  a(. . .)  are  related  to  the  Wigner  3-j  symbols  used  in  quantum 
mechanics. 

For  an  array,  the  translation  formula  can  be  used  in  conjunction  with  a  “T-matrix”  (T  for 
transition)  method  to  determine  the  total  pressure  field  of  an  array  of  scatterers/radiators 
relative  to  theoretically  any  coordinate  system.  It  is  especially  useful  in  cases  for  which 
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multiple-scattering  and  near  field  pressure  calculations  are  important.  A  good  general  refer¬ 
ence  on  the  subject  is  the  proceedings  of  an  ONR-sponsored  international  symposium  held 
in  1979  [6].  The  method  was  originally  developed  to  solve  acoustic  and  electromagnetic 
scattering  from  a  single  obstacle.  Later,  it  was  applied  to  the  scattering  of  elastic  waves, 
and  to  scattering  from  arrays  of  obstacles.  It  is  important  to  note  that  the  T-matrix  method 
accounts  for  multiple  scattering  to  all  orders,  and  that  the  method  is  capable  of  solving  array 
problems  with  arbitrarily  dense  -  randomly  packed  elements. 

The  scattering  properties  of  each  unique  scatterer  (transducer)  in  an  array  are  described 
by  a  so-called  T-matrix,  using  a  discrete  basis  set  of  spherical  harmonics.  The  T-matrix  for 
a  particular  finite  obstacle  is  determined  by  analyzing  the  scattering  characteristics  of  the 
obstacle  subject  to  “standing”  spherical  waves  of  the  form 

pl(r,e,<f>)=jn(kr)W(074>) 

The  form  of  these  incident  pressures  is  dictated  by  the  spherical  addition  formula,  and  the 
ability  to  represent  an  arbitrary  incident  pressure  field  in  terms  of  such  a  basis.  The  diagram 
below  graphically  indicates  the  T-matrix  use.  The  Rn  are  essentially  reflection  coefficients 
from  the  obstacle,  and  constitute  entries  into  the  T-matrix. 

T-matrix 

E  W*»0  W  <j>)  ^  E  Rnhn(kr)a%(6,  <f>) 

n  n 

N - - - '  ^ - s, - ' 

standing  waves  radiating  waves 

To  demonstrate  the  use  of  the  addition  formulas  with  the  T-matrix  formalism,  consider  a 
two  element  array  in  which  both  transducers  are  active.  The  two  translations  of  coordinates 
needed  can  be  found  simultaneously  due  to  the  simple  relationship  between  angles  of  one 
system  relative  to  the  other.  Let  G,j  represent  this  translation  formula  from  system  j  to 
system  i  (by  applying  the  spherical  addition  formula).  Furthermore,  let  Tk  be  the  “T- 
Matrix”  for  the  kth  transducer  modeled,  and  the  radiated  amplitudes  of  the  two  transducers 
(in  the  absence  of  their  neighbors)  be  given  by  the  vectors  Bx  and  B2.  One  obtains  the 
following  system  of  equations  for  the  unknown  scattering  pressure  amplitudes  A\  and  A2: 

TiGi2(A2  +  B2)  =  A\ 

T2G2i(Ax  +  B\)  =  A2 

or  in  matrix  form: 

(  I  -Ti  Gj2  \  f  A\  \  _  /  Tx  G12B2  \ 

V  -T2  G21  1  )\A 2  )  V  T2  G 21Bj  J 


(  0  Tx  G j  2  \ 

\  T2  G2i  0  J  ’ 


If  we  let 
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the  matrix  equation  becomes  even  simpler: 


(1-5) 


Once  the  scattering  amplitudes  are  found,  the  total  pressure  at  any  field  point  may  be 
found  by  first  representing  that  field  point  in  terms  of  each  individual  coordinate  system 
— >  ( rj,0j,4>j )  and  summing  the  series  over  each  array  element: 


N  n 

p(r,  0,  <j>)  =  J2  12  12  \Aimn  +  Bjmn]hn(kr j)£l™  (0j,  <f>j) 

j  n— 0  77i— — n 

A  particular  application  of  the  above  procedure  is  that  of  determining  source  levels  for 
an  array  of  16  closely  packed  spherical  shell  transducers.  Let  the  array  consist  of  two  rows  of 
8  spherical  shell  transducers  (with  outer  radius  “a”  equal  0.5  meters)  driven  at  a  frequency 
such  that  ka  =  1,  and  having  separation  distances  of  7t/2  =  A/2  «  1.57  meters  along  each 
row  and  7r/4  =  A/4  «  .785  meters  between  rows.  Furthermore,  one  row  of  transducers  is  90° 
out  of  phase  with  that  of  the  other.  A  graphic  of  the  array  is  given  below: 


wavelength 
h=K  m 


shell  radius  (a):0.5m 

iJll  shell 
{IgL  thickness:  .01  m 

/UpP  c  1=1 490  m/s 

5p  ,=474Hz 

t!  ka=l 


%/4m 

Orientation  of  the  spherical  shells 
in  the  16  element  array 


By  using  the  T-matrix  for  a  spherical  shell,  and  the  translation  formulas,  the  far  field 
source  level  can  be  found  for  the  array  operating  at  a  frequency  of  474  Hertz  which  corre¬ 
sponds  to  an  acoustic  wavelength  of  7t/2  (ka  =1).  The  figure  below  is  the  far-field  pattern 
for  the  array,  scaled  by  the  far  field  amplitude  of  a  single  spherical  shell  (or  source  level  of 
a  single  element  of  the  array).  The  graph  displays  source  levels  in  dB  in  the  plane  of  the 
array  as  well  as  a  superimposed  image  of  the  array  orientation.  Angles  are  measured  from 
the  positive  z  axis. 
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Source  level  of  array,  ka=1 ,  analytic  scaled  by  FF  of  a  single  spherical  shell 


270 


2  Validation  Procedure 

To  test  the  effectiveness  of  a  boundary  element  or  a  finite  element  code  to  construct  the  T- 
matrix  of  a  particular  transducer,  comparison  is  made  to  an  idealized  physical  problem  for 
which  there  is  an  analytical  solution.  The  idealized  problem  addressed  is  that  of  a  fluid  loaded 
thin  spherical  shell.  This  is  a  canonical  problem  which  lends  itself  to  analytic  treatment  due 
to  the  simple  characterization  of  the  boundaries,  and  the  fact  that  fluid  loading  does  not 
affect  mode  shape  (it  does  affect  the  eigenfrequencies  of  the  shell  however).  The  analytical 
solution  can  be  represented  in  terms  of  spherical  harmonics,  and  the  resulting  T-matrix  is 
diagonal.  Following  Junger  and  Feit  [2]  the  modal  mechanical  impedance  of  a  thin  spherical 
shell  is: 
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a  shell  radius 

h  shell  thickness 

E  Young’s  modulus 

v  Poisson’s  ratio 

Vt  ua/cp 

cp  plate  wave  speed  (\jE/ps(l  —  *)) 

cj  fluid  wave  speed 

u  frequency 

ps  ,  p/  solid  and  fluid  densities 
n  index  of  the  spherical  harmonic 

(It  should  be  noted  that  the  above  impedances  are  independent  of  the  index  “m”  which 
characterizes  the  longitudinal  variation  of  the  field  quantities.)  From  the  definition  of  the 
modal  mechanical  impedance  we  have 

Inmjn (&Gt)  Rnmh"n{kE)  =  lLoZnWnm 

where  Inm  and  Rnm  are  incident  and  reflected  modal  amplitudes,  and  Wnm  is  the  modal 
normal  displacement  of  the  shell.  The  vector  form  of  the  amplitudes  requires  an  ordering 
of  the  unknowns,  which  is  as  follows:  fi*o,o  Ri  ,  Ri,-i  -4  R2  ,  Ri,o  ->•  R3  ,  R\,\  — ► 
R4  ,  R-2,-2  R5  ,  R'2,-1  Re  ,  —  A  second  relationship  between  the  pressures  and 
displacements  at  the  surface  of  the  shell  results  from  Euler’s  equation  (the  assumption  of  no 
cavitation  on  the  surface  of  the  shell): 

Inmjn  (A:fl)  T  Rnmh"n  (fco)  — 

Combining  these  and  eliminating  the  normal  displacement  one  can  solve  for  the  scattered 
amplitude  in  terms  of  the  incident  amplitude: 

iZjJr/jka)  ~  pjCfjn(ka)  ' 

— iZnhn'(ka)  +pfCfhn(ka)_  nm 

The  bracketed  term  above  is  an  analytical  expression  for  the  diagonal  entries  of  the  T-matrix 
for  a  thin  spherical  shell. 

2.1  ATILA  modeling  for  T-matrix  determination 

T-matrix  entries  have  been  tabulated  for  a  spherical  shell  using  ATILA  (version  5.11)  with 
the  following  set  of  input  parameters  (these  are  the  same  values  as  those  used  in  the  16 
element  array  given  in  the  preceding  section): 
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a 

h 

E 

v 

f  =  o;/2tt 
Ps  5  P/ 

CJ 


0.5  meters 
0.01  meters 
215  GPa 
0.33 

474  Hertz  ( ka  —  1) 

7500  Kg/m3  &:  1000  Kg/m3 
1490  meters/sec 


The  values  of  the  diagonal  entries  of  the  T-matrix  using  ATILA  are  given  in  Table  II.  The 
analytical  values  for  the  diagonal  entries  of  the  T-matrix  for  the  same  set  of  parameters  are 
given  in  Table  I. 


Table  I  -  T-matrix  diagonals  -  Analytic 

n 

Real  part 

Imag  part 

Ampl. 

Phase(degrees) 

i 

-0.13465E-01 

0.11525E+00 

0.11604E+00 

0.96663E+02 

2 

-0.60255E-02 

0.77390E-01 

0.77624E-01 

0.94452E+02 

3 

-0.60255E-02 

0.77390E-01 

0.77624E-01 

0.94452E+02 

4 

-0.60255E-02 

0.77390E-01 

0.77624E-01 

0.94452E+02 

5 

-0.61773E-03 

-0.24847E-01 

0.24854E-01 

-0.91424E+02 

6 

-0.61773E-03 

-0.24847E-01 

0.24854E-01 

-0.91424E+02 

7 

-0.61773E-03 

-0.24847E-01 

0.24854E-01 

-0.91424E+02 

8 

-0.61773E-03 

-0.24847E-01 

0.24854E-01 

-0.91424E+02 

9 

-0.61773E-03 

-0.24847E-01 

0.24854E-01 

-0.91424E+02 

Table  II  -  T-matrix  diagonals  .-  ATILA 

n 

Real  part 

Imag  part 

Ampl. 

Phase  (degrees) 

1 

-1.3397e-02 

1.3587e-01 

1.3652e-01 

9.5631e+01 

2 

-6.0482e-03 

8.1155e-02 

8.1380e-02 

9.4262e+01 

3 

-3.1 121e-03 

8.2552e-02 

8.261  le-02 

9.2159e+01 

4 

-6.0480e-03 

8.1158e-02 

8.1383e-02 

9.4262e+01 

5 

5.2933e-04 

-2.3758e-02 

2.3764e-02 

-8.8724e+01 

6 

8.0858e-04 

-2.3813e-02 

2.3827e-02 

-8.8055e+01 

7 

1.3738e-03 

-2.4142e-02 

2.4181e-02 

-8.6743e+01 

8 

8.0907e-04 

-2.3814e-02 

2.3828e-02 

-8.8054e+01 

9 

5.2948e-04 

-2.3759e-02 

2.3764e-02 

-8.8723e+01 

As  seen  in  the  tables,  the  ATILA  generated  elements  suffer  from  a  relatively  small  error 
in  the  phase  (within  3  degrees  of  the  analytic),  but  for  the  fundamental  mode,  the  amplitude 
values  are  off  by  about  18%  of  the  analytic  value  (for  the  remaining  terms  the  relative  error 
is  within  6%).  In  addition,  when  considering  the  off  diagonal  entries  of  the  ATILA  generated 
T-matrix,  there  appeared  to  be  “leakage”  of  scattered  energy  between  the  n  =  2,  m  =  2  and 
the  n  =  2,  m  =  —2  modes. 
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The  ATILA  run  made  use  of  the  “new”  curved  shell  elements.  The  spherical  shell  was 
divided  into  72  elements,  while  the  surrounding  fluid  contained  two  fluid  layers  which  were 
0.25  meters  thick  followed  by  eight  fluid  layers  one-half  meter  thick  for  a  total  fluid  thickness 
of  5  meters.  Each  fluid  layer  had  72  fluid  “block  elements”,  and  the  entire  finite  element 
mesh  was  truncated  with  dipolar  damping  elements  which  are  consistent  with  applying  the 
time-harmonic  Sommerfeld  radiation  condition: 

ikp(r,  0,<f>)  +  —  (r,  6,<f>)  +  -p(r,0, 4>)  — y  0  as  r oo 
or  r 

One  likely  source  of  error  is  due  to  the  approximate  nature  of  the  radiation  boundary  con¬ 
dition  applied  at  the  truncated  fluid  domain.  A  second  source  is  the  coarseness  of  the  mesh, 
which  has  at  the  poles  four  elements  with  longitudinal  separation  angles  of  ninety  degrees. 

The  amplitudes  of  the  radiating  harmonics  are  extracted  from  the  ATILA  generated 
scattered  pressure  at  given  field  points  by  application  of  a  second  code  which  employs  a 
singular  value  decomposition  of  the  total  scattered  field  into  modal  amplitudes.  In  these  runs, 
the  match  between  ATILA  scattered  pressure,  and  the  modal  representation  of  the  pressure 
is  extremely  close,  which  is  considered  in  checking  the  fidelity  of  the  modal  expansion  to  the 
raw  data. 

2.2  SYSNOISE  modeling  for  T-matrix  determination 

The  currently  available  version  of  SYSNOISE  is  release  5.3A.  LMS  technologies  has  agreed 
to  make  the  Naval  Postgraduate  School  a  beta  site  for  the  testing  of  version  5.4,  which  is 
anticipated  to  become  commercially  available  by  December  of  1998.  However,  as  of  October, 
we  have  yet  to  receive  a  copy  of  the  5.4  version.  This  version  is  necessary  since  it  allows  one  to 
use  input  from  ATILA  generated  impedance  matrices  for  the  determination  of  scattering  and 
radiation  from  piezoelectric  transducers.  SYSNOISE  version  5.3A  is  incapable  of  accepting 
inputs  from  ATILA.  For  these  reasons,  the  focus  has  been  on  benchmarking  the  current 
version  of  SYSNOISE  to  the  same  problem  given  above  -  that  of  the  scattering  characteristics 
of  a  submerged  spherical  shell. 

Documentation  of  SYSNOISE  5.3A  indicates  that  it  is  possible  to  have  user  defined  input 
pressures  through  a  special  “user”  library.  A  FORTRAN  code  has  been  written  to  provide 
incident  pressures  of  the  form  necessary  to  produce  the  requisite  T-matrix  elements  (standing 
spherical  harmonics  of  the  form  p‘(r,0,<£)  =  jn(kr)Q%(0, <f>)).  Unfortunately,  the  makefile 
necessary  to  recompile  the  SYSNOISE  source  code  is  incompatible  with  the  operating  system 
(HPUX10.2).  The  design  team  in  Belgium  has  been  in  contact  with  Dr.  Scandrett  to  try 
and  produce  a  viable  makefile  which  would  allow  compilation  and  linking  of  the  SYSNOISE 
source  code  with  user  defined  functions,  but  as  of  now,  have  been  unable  to  do  so. 

As  an  alternative  to  using  user  defined  incident  pressures  of  the  form  desired,  one  may 
use  a  series  of  incident  plane  waves.  This  requires  that  the  incident  plane  waves  be  decom¬ 
posed  into  spherical  harmonics,  which  can  then  be  used  in  combination  with  their  scattered 
spherical  harmonics  to  produce  the  T-matrix.  Such  a  methodology  can  be  written  in  matrix 
form  as  follows: 
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TB  =  S 

where  T  is  the  sought  after  T-matrix,  the  columns  of  B  are  the  amplitudes  of  the  incident 
plane  waves  decomposed  into  standing  harmonics,  and  S  are  the  decomposed  scattered  pres¬ 
sure  amplitudes  in  terms  of  radiating  spherical  harmonics.  The  equation  is  transposed  and 
solved  in  a  least  squares  sense  for  the  T-matrix  of  the  transducer. 

Unfortunately,  representation  of  an  arbitrary  incident  plane  wave  as  a  sum  of  standing 
spherical  harmonics,  may  take  several  hundred  terms,  while  the  scattering  amplitudes  from 
the  finite  obstacle  have  typically  only  a  few  non-negligible  components.  In  particular,  if  only 
the  first  N  harmonics  are  necessary  to  adequately  represent  the  scattered  pressure  from  the 
obstacle,  only  NxN  T-matrix  is  needed.  Because  of  the  necessity  of  properly  representing 
a  plane  wave,  we  typically  require  the  number  of  rows  in  the  matrix  B  to  be  large  (say  M 
where  M  »  N).  The  number  of  columns  in  each  of  the  matrices  B  and  S  depend  upon  the 
number  of  distinct  incident  plane  waves  needed  to  determine  the  T-matrix  (say  P  where  we 
must  have  P  >  N).  The  full  T-matrix  given  these  considerations  must  be  M  x  M,  but  we’re 
only  interested  in  the  upper  left  block  of  this  matrix  (7n).  We  have  the  matrix  equation 
rewritten  below 


Tu  Tn  \(B1\_(S1\ 

T21  T-n  ){b2  )-{  S2  ) 

where  the  blocks  T12,  T2\,  and  T22 ,  are  respectively  N  x  (M  —  N),  (M  —  N)  x  N,  and  (M  — 
N)  x  (M  —  /V),  respectively.  The  full  T-matrix  for  an  arbitrary  scatterer  could  theoretically 
be  found  using  only  plane  waves,  but  it  would  require  the  number  of  columns  in  S  and  B 
to  be  at  least  as  large  as  M  -  the  number  of  harmonics  needed  to  represent  the  incident 
plane  wave.  This  highlights  our  interest  in  having  the  ability  to  arbitrarily  specify  incident 
pressures  through  user  defined  functions. 

For  the  special  case  of  a  spherical  shell,  we  know  that  the  oflf  diagonal  entries  of  the 
T-matrix  should  be  zero,  and  hence  Tv2  =  T2\  —  0.  Exploiting  this  property,  we  can  rewrite 
the  first  N  rows  of  the  above  matrix  equation  in  the  form 

T\\B\  = 

Now  the  number  of  columns  in  the  S  matrix  need  only  be  greater  than  or  equal  to  N  the 
row/column  size  of  the  sought  after  T-matrix.  For  a  general  obstacle/transducer,  this  trick 
cannot  be  used  due  to  mode  coupling  of  the  spherical  harmonics. 

Before  the  above  methodology  is  exploited,  a  check  is  made  on  the  ability  of  SYSNOISE 
to  accurately  determine  the  plane  wave  scattering  from  a  fluid  loaded  structure.  To  this  end, 
the  plane  wave  scattering  from  the  same  spherical  shell  (using  the  same  input  parameters 
as  above)  is  found  and  compared  to  an  analytical  solution.  The  incident  plane  wave  in  this 
instance  is  propagating  in  the  positive  ‘z’  direction.  The  analytical  solution  can  be  written 
for  the  scattered  pressure  in  terms  of  radiating  spherical  harmonics: 
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ps(r,0,<t>)  =  Yj  Anhn(kr)Pn(cosO) 


n= 0 


An  =  (— i)”(2n  +  1)- 


- Znj'n(ka )  +  ipfCfjnika) 


Znh'n(ka)  -  ipjCjhn{ka) 

and  Zn  is  the  modal  mechanical  impedance  given  previously.  Comparison  of  SYSNOISE 
output  with  the  analytic  solution  is  shown  in  the  figures  below. 


Plane  wave  scattering,  SYSNOISE  &  analytic,  Freq=474Hz,  ka=1 


Z  coordinate  value 


Z  coordinate  value 

The  discretization  used  in  SYSNOISE  consisted  of  216  planar  quadrilateral  elements  to  model 
the  shell.  It  can  be  seen  that  the  relative  error  is  greatest  in  the  backscatter  direction,  and 
is  particularly  bad  if  field  points  too  close  to  the  shell  are  employed.  In  total  however,  the 
results  are  favorable,  being  by  and  large  within  a  relative  error  of  10  percent.  We  are  now 
in  a  position  to  try  and  determine  the  9x9  T-matrix  elements  of  the  spherical  shell  using 
9  incident  plane  waves(this  will  include  spherical  harmonics  through  the  quadrupole  term), 
and  their  resulting  scattered  pressures.  The  B  matrix  has  column  entries  (one  column  for 
each  plane  wave)  which  are  the  analytic  values  of  the  coefficients  used  in  expanding  the  given 
incident  plane  wave  in  terms  of  standing  harmonics  through  the  formula  [5] 

oo  n 

e-ilro“7  =  E  E 

n=0  m=  — n 
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Pn.m  — 


(— i)n(2n  +  l)(n  —  m)\ 


P™( cos  a)e 


—imj3 


(n  +  m)! 

C037  =  sin  a  sin  6  cos (<f>  —  0)  +  cos  a  cos  0 

and  the  values  of  a  and  0  are  the  azimuthal  and  longitudinal  angles  of  the  plane  wave’s 
propagation  direction  vector.  The  components  of  the  S  matrix  are  found  using  the  scattered 


Plane  wave  scattering,  SYSNOISE  &  analytic,  Freq=474Hz,  ka=1 


pressures  at  field  points  at  a  distance  of  5  meters  from  the  sphere  center  generated  by 
SYSNOISE,  and  analyzed  using  a  least  squares  program.  The  relative  residual  error  was  in 
every  case  on  the  order  of  6  x  10-5.  (This  relative  error,  is  the  length  of  the  error  vector  in 
comparing  the  SYSNOISE  generated  data  to  its  computed  spherical  harmonic  representation, 
divided  by  the  length  of  the  original  data  vector.  It  is  symbolically  represented  as 

,  .  \\S-TH\\ 

relative  error  = - = - 

Ill'll 

■*—4 

where  the  vector  H  corresponds  to  the  set  of  radiating  spherical  harmonics  at  each  of  the 
prescribed  field  points.) 


adjusted  FF  Pressure  phase  ,  relative  error  r  times  FF  Pressure  amplitude 
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Plane  wave  scattering,  SYSNOISE  &  analytic,  Freq=474Hz,  ka=1 ,  Far-field 


Plane  wave  scattering,  SYSNOISE  &  analytic,  Freq=474Hz,  ka=1 ,  Far-field 


The  final  results  of  the  T-matrix  determination  are  given  in  Table  III. 
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Table  III 

-  T-matrix  diagonals  -  SYSNOISE 

n 

Real  part 

Imag  part 

Ampl. 

Phase  (degrees) 

T 

-1.2453e-02 

1.1093e-01 

1.1163e-01 

9.6405e+01 

2 

-5.5062e-03 

7.3775e-02 

7.3981e-02 

9.4268e+01 

3 

-5.5539e-03 

7.3788e-02 

7.3996e-02 

9.4304e+01 

4 

-5.5532e-03 

7.3784e-02 

7.3993e-02 

9.4304e+01 

5 

-5.5229e-04 

-2.2295e-02 

2.2301e-02 

-9.1419e+01 

6 

-5.3767e-04 

-2.2251e-02 

2.2258e-02 

-9.1384e+01 

7 

-5.371 6e-04 

-2.2120e-02 

2.2127e-02 

-9.1 391e+01 

8 

-5.4100e-04 

-2.2240e-02 

2.2246e-02 

-9.1393e+01 

9 

-5.4614e-04 

-2.2257e-02 

2.2263e-02 

-9.1406e+01 

As  can  be  seen  in  the  above,  there  is  a  much  better  match  to  the  analytic  values  than 
was  found  in  using  ATILA  alone.  In  particular,  the  relative  errors  in  the  magnitude  are 
about  4  percent  for  the  breathing  mode,  nearly  5%  for  the  dipole  terms,  and  increase  to 
about  10  percent  for  the  quadrupole  modes.  The  phase  error  is  uniformly  less  than  one-half 
of  1  percent.  There  are  however,  off  diagonal  entries  in  the  computed  T-matrix  which  are 
non-negligible  in  the  above  analysis.  For  the  first  7  modes  of  excitation,  the  amplitudes 
of  T-matrix  entries  in  the  corresponding  column  are  roughly  two  orders  of  magnitude  less 
that  of  the  diagonal  entry.  For  the  n=2  m=l  and  m=2  columns,  it  is  found  that  the 
n=2,  m=-l  and  m=-2  have  amplitudes  only  one  order  of  magnitude  less  than  the  diagonal 
entry.  It  is  suspected  that  this  is  the  result  of  two  phenomenon.  The  first  of  these  is  that 
the  relationship  between  associated  Legendre  functions  tends  to  magnify  those  Legendre 
functions  with  negative  index  m  as  seen  by  the  identity 


p-m 


,  x  (n  —  m)!  „  ,  v 
w = (i) 


A  second  reason  for  the  larger  than  expected  magnitudes  of  these  off-diagonal  entries  is  that 
perhaps  these  modes  are  preferentially  magnified  by  the  mesh  used  (in  this  case  a  set  of  216 
flat  quadrilateral  shapes  generated  by  SYSNOISE). 


2.3  ATILA  modeling  with  EQI  for  T-matrix  determination 

In  an  entirely  analogous  fashion,  the  T-matrix  can  be  determined  from  ATILA  coupled  with 
the  EQI  code  which  is  a  boundary  integral  method  based  upon  Jones’  technique  [1]  for 
handling  anomalous  interior  eigenfrequencies.  The  discretization  used  for  the  full  ATILA 
code  was  employed  on  just  the  surface  of  the  spherical  shell,  and  coupled  to  EQI  to  find  the 
results  of  scattering  calculations.  This  work  was  accomplished  with  the  help  of  Regis  Bossut, 
(an  I  SEN  employee)  since  NPS  does  not  currently  have  the  rights  to  use  EQI. 

At  present,  EQI  can  not  handle  arbitrary  incident  pressures  in  its  boundary  integral 
formulation,  so  the  same  methodology  used  with  SYSNOISE  will  be  employed  here.  It 
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should  be  noted,  that  the  designers  of  ATILA/EQI  are  open  to  modifications  of  their  code 
which  would  allow  user  defined  incident  pressures  should  this  research  continue,  and  this 
report  supplies  adequate  justification  for  purchasing  such  a  development. 


ATILA  &  analytic, Freq=474Hz,  Ka=1 ,  PW  scatter  at  0.5  meter 


As  with  SYSNOISE,  the  EQI  boundary  element  code  is  tested  with  plane  wave  scattering 
from  a  spherical  shell.  The  first  graph  compares  the  analytic  scattered  pressure  to  the  EQI 
results.  The  asterisks  are  points  on  the  surface  of  the  scatterer  used  to  determine  field  point 
pressures  elsewhere  in  the  fluid.  As  can  be  seen,  these  are  quite  accurate,  but  in  using  these 
values  to  determine  the  field  point  pressures  on  the  surface  of  the  shell  would  lead  to  large 
errors.  A  second  graph  compares  the  analytic  solution  with  ATILA/EQI  using  field  points 
at  1  meter  from  the  center  of  the  shell  (1/2  meter  from  the  shell  surface).  At  this  distance, 
error  from  field  point  calculations  appears  to  have  been  avoided.  It  is  at  this  distance,  that 
field  point  scattered  pressures  will  be  computed  for  the  T-matrix  determination  given  the 
same  set  of  incident  plane  waves  used  in  the  SYSNOISE  analysis. 
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ATILA  &  analytic,  Freq=474Hz,  Ka=1,  PW  scatter  at  1  meter 


The  results  of  this  study  lead  to  the  following  table  of  diagonal  entries  for  the  T-matrix. 


Table  IV 

-  T-matrix  diagonals  -  ATILA/EQI 

n 

Real  part 

Imag  part 

Ampl. 

Phase(degrees) 

i 

-1.2869e-02 

1.1277e-01 

1.1350e-01 

9.6510e+01 

2 

-6.4242e-03 

7.9990e-02 

8.0248e-02 

9.4592e+01 

3 

-6.4859e-03 

8.0094e-02 

8.0356e-02 

9.4630e+01 

4 

-6.4356e-03 

7.9975e-02 

8.0234e-02 

9.4601e+01 

5 

-5.9398e-04 

-2.4592e-02 

2.4599e-02 

-9.1384e+01 

6 

-6.0821e-04 

-2.4650e-02 

2.4657e-02 

-9.1413e+01 

7 

-6.2987e-04 

-2.5082e-02 

2.5090e-02 

-9.1439e+01 

8 

-6.0554e-04 

-2.4649e-02 

2.4656e-02 

-9.1407e+01 

9 

-5.9399e-04 

-2.4592e-02 

2.4600e-02 

-9.1384e+01 

In  the  above  table  we  have  the  best  results  so  far.  The  monopole  term  has  a  relative  error 
of  about  2%,  the  dipole  about  4%,  and  the  quadrupole  about  1%.  It  is  believed  that  a  direct 
cause  of  the  larger  relative  error  in  the  dipole  term  is  due  to  the  rather  coarse  discretization 
at  the  two  poles  of  the  surface  mesh.  The  EQI  result  also  suffered  less  from  energy  leakage 
in  the  quadrupole  terms  than  SYSNOISE  or  ATILA  alone,  but  there  is  still  some  occurring 
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from  the  n=2,  m=2  mode  to  the  n=2,  m=-2  mode.  In  the  seventh  column,  the  entry  in  the 
fifth  row  (n=2,m=-2)  is  only  one  order  of  magnitude  less  than  the  diagonal  entry. 

2.4  Comparison  of  methods  in  Array  calculations 

Using  the  model  16  element  array,  for  which  we  have  an  analytically  determined  solution,  we 
compare  each  of  the  three  computed  T-matrices.  The  computed  9x9  T-matrices  are  read 
into  the  array  code,  and  a  comparison  of  source  levels  and  surface  pressures  were  undertaken. 
The  polar  plot  of  the  array  source  level  was  unable  to  distinguish  the  four  curves,  and  so 
is  not  repeated  here.  Instead,  line  graphs  indicated  the  amplitude  (absolute  pressure  and 
pressure  in  dB)  and  phase  of  the  different  solutions  are  given  below. 

For  the  far-field  calculations  in  dB,  we  see  that  the  three  methods  all  do  quite  well.  The 
error  in  dB  is  uniformly  less  than  2  for  all  methods,  with  ATILA/EQI  best,  and  ATILA  alone 
worst.  The  results  are  scaled  by  the  magnitude  of  a  single  shell  radiating  in  a  breathing  mode. 
It  can  be  seen  in  the  graph,  that  the  ATILA/EQI  results  are  within  1  dB  of  the  analytic, 
and  for  much  of  the  angles  where  the  dominant  lobes  of  the  pressure  amplitude  reside,  the 
error  is  well  below  1  dB. 


Array  far-field  pressure  level  in  dB,  ka=1 


An  analogous  graph  using  the  absolute  pressure  in  the  far  field  clearly  indicates  the 
differences  in  solutions.  Note  that  the  circles  in  the  lower  graph  represent  angles  from 
which  incident  plane  waves  were  coming  in  the  T-matrix  determination  for  SYSNOISE  and 
ATILA/EQI.  However,  except  for  the  180°  circle,  the  incident  plane  waves  were  not  propa- 
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gating  in  the  plane  of  the  array,  which  is  the  plane  for  which  the  results  have  been  reported. 

Source  level  of  array,  analytic  vs  sysnoise  vs  eqi  vs  atila 
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in  addition  to  the  far-field  source  level  values  of  a  given  array.  To  test  this  ability,  two  of  the 
spherical  shells  near  the  acoustic  center  of  the  array,  were  analyzed  for  their  surface  pressure 


Total  surface  pressures  in  plane  of  array  for  shell  #8 
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(amplitude  and  phase).  They  are  the  two  shells  numbered  8  and  9  in  the  array  graphic  given 
previously.  Spherical  shell  8  is  on  the  “shadow”  side  of  the  array  while  9  is  on  the  “lit”  side. 

Total  surface  pressures  in  plane  of  array  for  shell  #9 
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The  graphs  indicate  an  error  (relative  to  the  max  pressure  amplitude)  which  is  less  than 
5%  for  nearly  the  whole  surface.  The  analytic  values  of  the  total  surface  pressure  used 
spherical  harmonics  up  to  a  value  of  n=6  for  the  array.  Only  the  total  surface  pressure 
in  the  plane  of  the  array  (y-z  plane)  are  shown,  but  this  is  the  plane  where  the  anticipated 
largest  and  smallest  values  of  the  total  pressure  should  occur.  In  fact,  the  pressure  max  points 
appear  to  lie  on  the  north  and  south  poles  of  the  shells  as  one  would  expect,  since  along 
these  directions  there  are  more  shells.  An  interesting  result  is  the  surface  pressure  phase  on 
shell  9.  In  this  instance,  ATILA/EQI  appears  to  do  the  worst  of  all  three  methodologies. 
It  is  assumed  that  this  result  is  due  to  the  relative  error  in  calculating  the  dipole  T-matrix 
elements 

The  important  thing  to  note  is  the  accuracy  of  the  actual  surface  pressure  for  this  array 
which  has  shells  of  radius  50  cm,  and  nearest  neighbors  only  11  cm  away. 
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